package com.record.utils;

import org.apache.commons.math3.distribution.NormalDistribution;
import java.util.Arrays;

/**
 * Shapiro–Wilk 正态性检验
 * 输出 W 和 p-value
 * 要求：样本量 ≥ 3
 */
public class ShapiroWilkTest {

    private double W;
    private double pValue;

    public ShapiroWilkTest(double[] data) {
        calculate(data);
    }

    public double getW() {
        return W;
    }

    public double getPValue() {
        return pValue;
    }

    /**
     * 主计算逻辑
     */
    private void calculate(double[] data) {
        int n = data.length;
        if (n < 3) {
            throw new IllegalArgumentException("样本量必须 ≥ 3");
        }

        // 排序
        double[] x = Arrays.copyOf(data, n);
        Arrays.sort(x);

        // 均值
        double mean = Arrays.stream(x).average().orElse(0);
        // 总平方和
        double ssd = 0.0;
        for (double v : x) {
            ssd += Math.pow(v - mean, 2);
        }

        // ✅ 如果所有值相同，方差为 0，无法计算，直接返回 W=1, p=1
        if (ssd == 0) {
            W = 1.0;
            pValue = 1.0;
            return;
        }

        // 理论正态分位数
        double[] m = new double[n];
        NormalDistribution nd = new NormalDistribution(0, 1);
        for (int i = 0; i < n; i++) {
            m[i] = nd.inverseCumulativeProbability((i + 1.0 - 0.375) / (n + 0.25));
        }

        // 标准化 m
        double mm = Arrays.stream(m).map(v -> v * v).sum();
        for (int i = 0; i < n; i++) {
            m[i] /= Math.sqrt(mm);
        }

        // 计算分子
        double numerator = 0.0;
        for (int i = 0; i < n / 2; i++) {
            numerator += m[i] * (x[n - 1 - i] - x[i]);
        }

        // 计算 W
        W = (numerator * numerator) / ssd;

        // 计算 p 值（使用近似公式）
        pValue = approximatePValue(W, n);
    }

    /**
     * p-value 近似估算函数
     */
    private double approximatePValue(double W, int n) {
        // 避免出现 log(负数)
        if (W <= 0 || W >= 1) {
            return (W >= 1) ? 1.0 : 0.0;
        }

        double y = Math.log(1 - W);
        // 根据 Royston (1995) 近似
        double mu = -1.2725 + 1.0521 * (Math.log(n) - 3.0);
        double sigma = 1.0308 - 0.26758 * (Math.log(n) - 3.0);
        double z = (y - mu) / sigma;

        NormalDistribution nd = new NormalDistribution(0, 1);
        return 1 - nd.cumulativeProbability(z);
    }

    /**
     * 测试入口（可独立运行）
     */
    public static void main(String[] args) {
        double[] data1 = {1, 1, 1, 1};           // 全部相同
        double[] data2 = {1, 1.1, 0.9, 1.05, 0.95}; // 近似正态

        ShapiroWilkTest test1 = new ShapiroWilkTest(data1);
        ShapiroWilkTest test2 = new ShapiroWilkTest(data2);

        System.out.println("【样本1】");
        System.out.println("W = " + test1.getW());
        System.out.println("p-value = " + test1.getPValue());
        System.out.println();

        System.out.println("【样本2】");
        System.out.println("W = " + test2.getW());
        System.out.println("p-value = " + test2.getPValue());
    }
}
